% Copyright (c) 2012, John Gebbie. All rights reserved.
%
% This library is free software; you can redistribute it and/or modify it
% under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation; either version 2.1 of the License, or (at
% your option) any later version.
%
% This library is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser
% General Public License for more details.
%
% You should have received a copy of the GNU Lesser General Public License
% along with this library; if not, write to the Free Software Foundation,
% Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301  USA


% Script for testing varying bathymetry.


%% Init

clear;
clear classes;
clc;

obj = ortrace();

%% Parameters

max_range = 400;

params = [];

% SSP
iRng = 1;
params.ssp_cz{iRng}.c = [ 1500 1510 ];
params.ssp_cz{iRng}.z = [ 0 500 ];
params.ssp_r(iRng) = 0;
params.ssp_step_z = 10;

% SEABED
params.seabed_range = 0:1:max_range;
params.seabed_depth = 300+200*sin(2*pi/200*params.seabed_range);
params.seabed_c = repmat(1600, size(params.seabed_range));
params.seabed_rho = repmat(1.6, size(params.seabed_range));

% RAY
params.rs = 0;
params.zs = 5;
% ray angles (positive angles points toward seabed)
params.theta_mn_mx_n = [ [1 1]*27.5*pi/180 1 ];
params.output_ranges_mn_mx_n = [10 max_range 100];
% params.output_ranges_mn_mx_n = [100 100 1];
params.output_full_ray = false;
params.output_grazing_angle = true;
params.step_m = 1;
params.max_db_loss = 150;
params.max_depth = Inf;

%% Configure

tic
obj.init_environment(params);
toc

%% Cast Ray

tic
% profile clear; profile on;
obj.cast_ray_fan();
% profile off; profile viewer;
toc

%% Plot

ray_out = obj.ray_out;

clf;

% light gray lines - output ranges
tmp_x = linspace(params.output_ranges_mn_mx_n(1),...
    params.output_ranges_mn_mx_n(2),...
    params.output_ranges_mn_mx_n(3));
for ix = 1:length(tmp_x)
    xx = [1 1]*tmp_x(ix);
    yy = [max(params.seabed_depth) 0];
    hPlot = plot(xx,yy,'Color',[1,1,1]*.8);
    hold on;
end
hold on;

% seabed
plot(params.seabed_range, params.seabed_depth, 'k');

% surface
plot([0, params.output_ranges_mn_mx_n(2)], [0 0], 'k');

% output grid points
if numel(ray_out) > 0
    xx = ray_out(ortrace.POS_R, :);
    yy = ray_out(ortrace.POS_Z, :);
    plot(xx, yy, '-b');
end

% plot grazing angles for bottom bounce
xx = ray_out(ortrace.POS_R,:);
yy = ray_out(ortrace.POS_GANG,:)*180/pi;
mm = ray_out(ortrace.POS_BNCB,:) ~= 0;
xx = xx(mm);
yy = yy(mm);
plot(xx,yy,'ro');

% config
set(gca, 'YDir', 'reverse');
xlim([0 params.output_ranges_mn_mx_n(2)]);
title('ray');
set(gca, 'YTick', 0:30:max(params.seabed_depth));
axis equal;
axis tight;
